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Abstract 

This paper summarises a number of new, potentially significant, results, obtained 
recently by the author and his collaborators, which impact on various issues related 
to the gravitational iV-body problem, both Newtonianly and in the context of general 
relativity. 



1 Introduction and motivation 



The overall objective of the research reported herein is the application of ideas and 
techniques from modern nonlinear dynamics and nonequilibrium statistical mechanics 
to self-gravitating systems, both Newtonianly and in the context of general relativity, 
with particular emphasis on the gravitational iV-body problem. The basic motivation 
for this research is a desire to identify some of the physical processes which can play 
a role in determining the structure and evolution of self-gravitating systems. The 
results described here will, for specificity, typically be formulated in the language 
of galactic dynamics. However, it should be evident that they also have potential 
implications in other settings as well, including, e.g., statistical quantum field theory 
in the early Universe. 

It should be stressed that the principal focus here is not on the detailed modeling 
of any specific class of astronomical object where, in particular, other nongravitational 
effects, such as dissipative hydrodynamics, can be important. However, the results 
reported here should find relatively direct applications to the study of systems like 
elliptical and lenticular galaxies, which are known to be gas poor, albeit not as gas 
poor as they were ten years ago. 

In approaching the study of self-gravitating systems, there are several different 
approaches which one might adopt. At the most fundamental level, one can attack 
the full iV-body problem, either by performing and analysing numerical simulations 
or by proving (hopefully useful) mathematical theorems which provide insights into 
the qualitative character of the evolution. In either case, the focus here is not on solar 
system type problems, involving a relatively small number of masses, one or two of 
which are much larger than the others. Rather, the principal focus is on collections 
of large numbers N of objects, comparable in mass, in particular the problem of the 
iV — > oo limit. 

Conventional wisdom says that, in this large N limit, such a collection of compa- 
rable masses can be described by a collisionless Boltzmann equation, i.e., what the 
mathematician would term the Vlasov-Poisson or Vlasov-Einstein system. Such a 
description involves the assumption that, Newtonianly, particles follow trajectories in 
the self-consistent gravitational potential associated with the average mass distribu- 
tion. Relativistically, one supposes that the particles follow geodesies in a spacetime, 
the form of which is determined by the stress energy tensor associated with the av- 
erage mass distribution. In this general context, two things need to be done, namely 
(1) to determine precisely the conditions under which such a mean field description 
is justified and then (2) to understand the qualitative and quantitative implications 
of this description. 
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In this connection, it is important to stress that, although the Vlasov-Poisson 
system was first formulated by Jeans in the context of galactic dynamics Jj]] more than 
twenty years before it was formulated in plasma physics, [0] the gravitational equation 
is substantially less well motivated than is the plasma analogue. In particular, the 
absence of shielding prevents a systematic derivation except for the special case of 
cosmological systems, which are nearly homogeneous and of finite age.|3| The problem 
becomes especially acute for the case of relativistic systems, the point being that the 
derivation of the corresponding Vlasov-Maxwell system relies heavily on the linearity 
of Maxwell's equations,^] whereas Einstein's equation is nonlinear. 

A yet simpler tact involves the consideration of orbits in a fixed potential. The 
idea here is to specify one's favourite potential, chosen to represent the bulk potential 
of some physical object, to study in detail the form of orbits in this potential, and only 
at the end of the day to incorporate the fact that the potential must be determined 
self-consistently by the mass distribution associated with the orbits themselves. 

If one chooses to focus simply on an average potential, in the context of a Vlasov 
description or otherwise, there remains the important question of determining pre- 
cisely, both qualitatively and quantitatively, what sorts of effects have been ignored. 
In other words, there remains the problem of structural stability. Only to the extent 
that these additional effects have been appropriately identified and understood can 
one say that one has a satisfactory understanding of so-called "collisionless stellar 
dynamics." 

The aim of this talk is to illustrate each of the preceding aspects by explaining 
several new results that have been derived during the past three or four years. Sec- 
tion 2 focuses on the full Newtonian iV-body problem. Section 3 then turns to the 
collisionless Boltzmann equation. Section 4 addresses several issues related to the 
problem of orbits in a fixed potential, and Section 5 concludes with a discussion of 
some aspects of the problem of structural stability. 

2 The gravitational iV-body problem 

Detailed numerical simulations over the past thirty years have established that, given 
generic initial conditions corresponding to a bound configuration, a self-gravita-ting 
system of point masses will typically exhibit a rapid approach, on a characteristic 
crossing time t cr , towards a statistical quasi-equilibrium where, in a time-averaged 
sense, the system only shows subsequent systematic variability on substantially longer 
timescales. Moreover, there has evolved a substantial and detailed conventional wis- 
dom which serves, at least roughly, to determine what kinds of initial data give rise to 
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what kinds of final quasi-equilibria. However, despite these impressive successes, one 
does not completely understand this process. Indeed, at a truly fundamental level 
there is no real understanding of why this happens. 

To obtain some insights into this basic question, it is natural to identify various 
microscopic and mesoscopic phenomena which could perhaps conspire to produce the 
qualitative macroscopic evolution that is observed in numerical experiments. The aim 
here is to identify two such phenomena. 

1. Viewed microscopically in the many particle phase space, Newtonian N-body sim- 
ulations exhibit an exponentially sensitive dependence on the specific choice of initial 
conditions. Specify unperturbed initial data, {r^(0) p^(0)}, (A = 1, ...,N), and per- 
form a simulation. Then specify perturbed initial data, {r^(0) p^(0)}, and repeat. 
What one discovers thereby is that, generically, quantities like the total iV-particle 
configuration space perturbation 

N N 

£l^(t)| 2 = £|r^)-r^)| 2 (1) 

A=l A=l 

will typically grow exponentially in time t on a relatively short time scale, even though 
the unperturbed and perturbed evolution are essentially identical at the macroscopic 
level. 

This fact was first observed by S. Ulam in the early 1960's, and the classic paper 
on the subject is by Miller. || However, only in the last several years, with the advent 
of improved computers, has this instability been studied systematically in complete 
and gory detail. |, g |, |, 

The net result of these investigations is that this instability is an exceedingly ro- 
bust phenomenon, which proceeds generically on a characteristic crossing time t cr , 
independent of many/most details. In particular, the timescale associated with this 
instability is independent of the detailed choice of initial data and the detailed choice 
of the initial perturbations, as well as the specific diagnostics used to quantify the 
growth of the perturbation - e.g., configuration or momentum space perturbation, 
the total iV-particle perturbation or the perturbation of "typical" or "representa- 
tive" particles. The timescale is also insensitive to a possible distribution of masses, 
provided that everything is not dominated by a few particularly massive particles. 

More significantly, the simulations also suggest strongly that the rate is insensi- 
tive to the total particle number N, provided at least that N ^ 2. Thus, e.g., for 
200 < N < 4000 one observes no appreciable changes if everything is scaled in terms 
of an iV-dependent characteristic time t cr . In particular, there is no sense in which 
the instability appears to "turn off" in the limit of large N. Indeed, Goodman et a/|9| 
have predicted that the instability should accelerate for large N, with the character- 
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istic timescale i* scaling as t* ~ t cr /lnN. Interestingly, t cr is the same timescale on 
which generic initial data evolve towards a macroscopic quasi-equilibrium. 
2. Gravitational N -body simulations evidence a considerable "memory." Viewed mi- 
croscopically or mesoscopically, the quasi-equilibrium does not correspond to a par- 
ticularly well-shuffled state. Many aspects of the initial data for individual particles 
are forgotten, but other aspects are in fact remembered. The strongest correlation 
between initial and final conditions, from which all others may perhaps derive, is 
between the initial and final values of the binding energy. Both for isolated sys- 
tems approaching a quasi-equilibrium and for collisions between pairs of objects, e.g., 
spherical polytropes and other axisymmetric or triaxial near-equilibria, there is a 
strong correlation between the initial and final values of the binding energy. 

Specifically, one observes that particles initially with high binding energy tend to 
end up with high binding energy, low with low, and intermediate with intermediate, 
even for an evolution that involves rapid, violent changes in the bulk potential, so that 
there is no obvious sense in which the binding energy should behave as an adiabatic 



invariant. [11, 12, W% 14 1 This phenomenon can be quantified at a coarse-grained 



mesoscopic level, through various binnings of the orbital data. [T^, Q Alternatively, 
it can be quantified at the completely microscopic level through the computation 
of a rank correlation between the initial and final orderings of particles in terms of 
their binding energies. [fl4| Such a computation shows that there exist strong, albeit 
not complete, correlations between initial and final conditions. The absence of a 
complete correlation is at least partly "collisional" in origin, but appears to persist 
even in the limit of large N, where, according to conventional wisdom, the system 
should be essentially collisionless in character. 

To summarise: In the gravitational iV-body problem one is confronted with a sys- 
tem that exhibits a rapid macroscopic evolution towards a statistical quasi-equilibrium. 
Viewed microscopically, this evolution is characterised by an exponentially sensitive 
dependence on the specific choice of initial conditions. However, despite this sensi- 
tive dependence, the quasi-equilibrium does not correspond, either microscopically or 
mesoscopically, to a particularly well-shuffled state. 



3 The collisionless Boltzmann equation 

The principal message of this section is that the collisionless Boltzmann equation, i.e., 
a mean field Vlasov description, is a constrained Hamiltonian system. This fact was 
first established by Morrison ||15|| for the electrostatic Vlasov-Poisson system, and is 
equally true for the analogous gravitational Vlasov-Poisson system. In this setting, the 
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fundamental dynamical variable is the one-particle distribution function, /(x, p, m), 
evaluated at a fixed time t, which is interpreted as a number density of particles of 
mass m at the spatial point x with conjugate momentum p. The gravitational po- 
tential $(t) is viewed as a functional of /(£), constructed in terms of an appropriate 
Green function. The phase space is then the infinite-dimensional space of distribu- 
tion functions. The Hamiltonian character of the evolution is manifest through the 
identification of a Hamiltonian function H[f] and a cosymplectic structure { . , . }, 
in terms of which the Vlasov equation takes the form df/dt = {/, H}. 

A generalisation to the spherically symmetric Vlasov-Einstein system is completely 



straightforward. ||16| The crucial point is that, for spherical systems, the gravitational 
degrees of freedom are not triggered: given appropriate boundary conditions, one can 
view the spacetime metric g ab at any given time t as a functional of the distribution 
function / at that same t. 

The full Vlasov-Einstein system is substantially more complicated, since, in the 
general case, one must incorporate the field degrees of freedom. However, the analysis 
still turns out to be straightforward, at least in principle. [17] The basic formulation 



is analogous to the Hamiltonian formulation of the Vlasov-Maxwell system, |18|, [19 
the nonlinearity of the Einstein equation not playing a significant role. 

Working in the context of the ADM formulation of general relativity, there are now 
three different dynamical variables, each defined on t = const hypersurfaces, namely 
(1) the distribution function, /(x, p, m), (2) the spatial three-metric, h ab (x), and 
(3) the conjugate field momentum, n afc (x). The natural phase space is the infinite- 
dimensional space coordinatised by these three variables. 

In this case the cosymplectic structure given as the sum of two pieces, namely 
(1) the functional Poisson bracket of vacuum gravity and (2) the matter bracket 
appropriate for the spherically symmetric Vlasov-Einstein system (which coincides 
also with the bracket for the Vlasov- Poisson system). Explicitly, for two functions 
F[f,h ab ,U ab ] and G[f,h ab ,U a % 



5F 5G 
Jf'Jf 



(2) 



where S/SX denotes a functional derivative with respect to the variable X and 



' dx l dpi dx l dpi 

denotes an ordinary three-dimensional Poisson bracket. 

To give meaning to variations 5X, one requires a rule identifying particle coordi- 
nates {x a ,p a } and {x' a ,p' a } in two nearby cotangent bundles. In the context of this 
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3 + 1 formulation, it is natural to identify spatial coordinates and conjugate momenta, 
as well as time t and mass m, i.e., 



x' = x p' = p t! = t m! = m. (4) 



However, other choices are also possible. [gO], ^TJ 

The Hamiltonian H = Hq + Hm = J d^xTic + I d?xli,M is also given as the sum 
of two pieces, namely (1) the ADM Hamiltonian Hq of vacuum gravity, i.e., 

H G = ^Th^Ul -® R + h-\u ab u ab - -U 2 ) 



16n 
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~2N b [D a (h- 1 / 2 U ab )} + 2D a (h-^ 2 N b U ab )^ (5) 

and (2) a matter Hamiltonian Hm constructed as the integral of the local energy 
density, i.e., 

H M = JdTf£ = J Nh}l 2 d?xT\. (6) 

Here dT = d 3 xd 3 pdm denotes the covariant seven- dimensional volume element on a 
t = const hypersurface, D a a covariant derivative on the hypersurface, and £ (x, p, m) 
= \p t \ the particle energy. N and N a correspond respectively to the lapse function 
and shift vector. 

This formulation reproduces the Vlasov-Einstein system in the sense that the 
equation dF/dt = (F,7i) for arbitrary functions F[f, h ab , U ab ] is equivalent to the 
Vlasov-Einstein system in its usual 3 + 1 form: The distribution function / satisfies 
df/dt = [£, /], and dh ab /dt and dl\ ab /dt both satisfy the appropriate 3 + 1 Einstein 
equations sourced by T a b [f]. For the spherically symmetric case, with the metric 
viewed as a functional of /, the first term in the bracket (F, G) disappears and the 
Hamiltonian H reduces to the ADM mass, Hadm-, realised as a volume integral in 
the natural fashion facilitated by Schwarzschild coordinates. 

Such a Hamiltonian formulation of the collisionless Boltzmann equation is in fact 
an example of a much more general result. Specifically, consider any Hamiltonian the- 
ory for a system with multiple degrees of freedom, and then construct the associated 
mean field description appropriate in the limit that correlations amongst the degrees 
of freedom are negligible (i.e., a statistical description in which the full N- "particle" 
distribution function is approximated as factorising into a product of reduced one- 
"particle" distribution functions). There is then a precise mathematical sense in 
which the mean field description of this Hamiltonian system is itself Hamiltonian. p2| 



The fact that the collisionless Boltzmann equation, or any other mean field theory, 
is Hamiltonian is significant in that a Hamiltonian evolution is much more restricted 
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than a non-Hamiltonian evolution. However, of more pragmatic importance per- 
haps is the fact that this Hamiltonian formulation permits, for the first time, a clear 
geometric approach to the problem of stability for general equilibrium solutions to 
the Vlasov-Einstein system. Here an "equilibrium solution" {/ , h® b , n^ 6 } entails a 
stationary matter distribution, corresponding to a spacetime that admits a timelike 
Killing field. 

The crucial fact facilitating this geometric approach is that every such equilib- 
rium is an energy extremal, so that the first variation S^'H vanishes identically for 
perturbations {Sf, Shab, 5U ab } which satisfy the constraints. The field constraints 
are enforced by restricting attention to perturbed initial data for which 5H/5N = 
5H/5N a = 0. The matter constraints, a reflection of Liouville's Theorem, im- 
ply that the perturbed distribution function must be generated from the unper- 
turbed /o via a canonical transformation in terms of some generating function h, 
i.e., f + 6f = exp([/i, .])/„ = /(,+ [h, f ] + \ [h, [h, /„]] + 

The fact that the equilibrium is an energy extremal implies that stability hinges on 
the sign of the second variation 5^ 2 'H. Indeed, modulo infinite-dimensional technical- 
ities the situation is analogous to the problem of stability for mechanical Hamiltonian 
systems. If S^H > for all infinitesimal perturbations, linear stability is guaran- 
teed. Alternatively, if S^H is indeterminate, one cannot necessarily infer a linear 
instability, but one does expect at least nonlinear instability and/or instability in the 
presence of dissipation. f23|j 

Indeed, one can actually prove that generic rotating axisymmetric equilibria are 
always unstable towards dissipation, as provided, e.g., by the emission of gravitational 
radiation. This is the collisionless analogue of the theorem [f25|j that all rotating 
perfect fluid stars are unstable. Neither for collisionless nor collisional systems is 
there any guarantee that the timescale associated with this instability is sufficiently 
short to be of interest astronomically. However, it is significant that general relativity 
triggers a generic instability which, apparently, is insensitive to the form of the self- 
gravitating matter. The astronomical implications of this instability are currently 
under investigation. 

This general approach to stability can also be adapted to the consideration of 
steady-state equilibria, such as an homogeneous and isotropic Friedman cosmology, 
where it provides an interesting derivation of the Jeans instability. |26]] Viewed New- 
tonianly, such an expanding Universe corresponds in the "inertial" frame to a sys- 
tem characterised by a time-independent Hamiltonian which finds itself in a time- 
dependent steady state. This explicit time-dependence can be removed by effecting 
a time-dependent canonical transformation into the average "comoving" frame. This 
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transformation leads to a new time-dependent Hamiltonian H(t). However, the first 
variation S^H vanishes identically, and the second variation S^H can be shown to 
satisfy dS^H/dt < 0. A simple Liapounov argument therefore guarantees that the 
existence of negative energy perturbations S^H < for sufficiently long wavelengths 
must imply an instability: If the system be perturbed in such a fashion that 5^H < 0, 
the energy can only become more negative, so that the "distance" from equilibrium, 
as probed by the magnitude of S^H, can only increase. 

4 Transient Ensemble Dynamics 

In studying the properties of orbits in a fixed potential, it is natural to apply the 
technology of nonlinear dynamics, as has been developed over the past several decades. 
However, in many settings involving gravitational systems, the utilisation of this 
technology may require a new twist. The standard formulation of nonlinear dynamics 
typically involves a theory of asymptotic orbital dynamics, which focuses primarily, if 
not exclusively, on the long time behaviour of individual orbits. However, for many 
astronomical systems this may not be appropriate. 

For example, in terms of their natural timescale, galaxies are relatively young 
objects, only ~ 100 — 200 crossing times t cr in age, so that it is not at all obvious 
that an asymptotic t — > oo limit is well motivated physically. Thus, e.g., standard 
estimates of Liapounov exponents typically require integrations for times t > I0 4 t cr , a 
period that is orders of magnitude longer than the age of the Universe, t H . Moreover, 
it is arguably true that, in many astronomical systems, individual orbits are not the 
fundamental objects of interest. It is, for example, obvious that one cannot track 
individual orbits of stars within a galaxy. All that one can detect are properties 
like the overall brightness distribution which reflect the contributions of many stars. 
Similarly, it is evident that one must focus on collections of orbits if he or she wishes 
ultimately to address the problem of self-consistency. 

For these reasons, it would seem more natural to consider instead a theory of 
transient ensemble dynamics, which focuses on the statistical properties of ensembles 
of orbits, restricting attention exclusively to short timescales, t < tn, and recognising 
that much of the observed behaviour may be intrinsically transient in character. 

This distinction may not be all that important for integrable, or near-integrable 
potentials, which contain only regular orbits. However, there is no reason to assume 
that the bulk potential associated with a self-gravitating system is integrable, or even 
near-integrable, and there are indications from numerical experiments that objects 
like rotating elliptical galaxies and barred spiral galaxies may admit large numbers 
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of stochastic, or chaotic, orbits. [ 

Roughly, regular orbits correspond to orbits that have regular shapes and are 
characterised by simple topologies, e.g., box orbits in two dimensions with trajectories 
that resemble a Lissajous figure or tube orbits in three dimensions that are restricted 
to a region with the topology of a torus. By contrast, stochastic orbits are manifestly 
irregular in shape and appear to "run all over" phase space. Unlike regular orbits, 
stochastic orbits exhibit an exponentially sensitive dependence on initial conditions, 
as manifest by the fact that such orbits are characterised by a positive Liapounov 
exponent 

The crucial point to be illustrated below is that, at least when considering stochas- 
tic orbits, the transient ensemble perspective can prove extremely important. 

Consider, for example, the scattering of photons incident on a multi-black hole 
system; or similarly, consider a star moving in a nonspherical potential, supposing that 
the star is unbound but that, owing to the shape of the potential, it can only escape 
in certain directions. For each of these problems, one discovers that, in certain phase 
space regions, the direction and time of escape to infinity exhibits a very sensitive 
dependence on initial conditions, in fact a fractal dependence, this being an example 
of what the nonlinear dynamicist would call chaotic scattering. PD| 

Naively, one might anticipate that this complex microscopic behaviour would lead 
to an equally complex description when considering the evolution of ensembles of 
orbits. However, this is not always the case. At least in certain cases, one observes 



instead striking regularities, which lead to a simple scaling behaviour, |3]J that may 
actually be universal. []32| The very fact that the microscopic evolution is complex 
appears to be responsible for the fact that the macroscopic evolution is very simple. 

As a concrete example, consider orbits in the two-dimensional potential V(x, y) = 
\{x 2 + y 2 ) — ex 2 y 2 , holding fixed the value of the energy E and studying as a function 
of e the evolution of orbits initially localised in a small phase space region. For e 
below a critical value e\ = 1/(AE), escape is impossible energetically. For values of e 
slightly above ei, escape is possible energetically, but only very few particles escape 
on short timescales and the time of escape exhibits no striking regularities, except 
that the escape probability eventually appears to decay towards zero. However, for e 
above another critical value €2 (only determined numerically), one sees the onset of 
striking scaling behaviour: 

1) After the decay of initial transients, the escape probability per unit time approaches 
a constant value P DO (e), which is independent of initial conditions, i.e., the location 
of the phase space region from which the initial ensemble was chosen. Moreover, this 
rate scales as -Poo(e) ~ (e — €2)" for a critical exponent a. 
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2) For fixed size of the initial phase space region probed by the ensemble, the time 

required to converge to also depends on e and satisfies Too(e) ~ (e — e 2 )~^. 

3) For fixed e, the convergence time depends on the linear size R of the phase 
space region that was probed initially, satisfying T^R) ~ R~ s 

Moreover, to within statistical errors a — (3 — <5 ~ 0. In a certain sense, the qualitative 
change in behaviour at e = e 2 is like a phase transition, complete with a critical 
"slowing down" as the "order parameter" e — e 2 — > 0. 

As another example in which the transient ensemble perspective is important, 
consider the behaviour of ensembles of orbits of fixed energy E, evolving in a two- 
dimensional time-independent potential V(x,y), which admits both regular and sto- 
chastic orbits. If V(x, y) is bounded from below and diverges at infinity, the constant 
energy hypersurfaces will be compact, so that the notion of "equilibrium" is well 
defined. One might therefore anticipate that generic ensembles of initial conditions 
will evolve towards an invariant distribution corresponding to a statistical equilibrium. 

To test this hypothesis, one can select localised ensembles of initial conditions 
of fixed E, corresponding to stochastic orbits initially located far from any regular 
phase space regions, and then evolve these initial data into the future. At least for 
certain potentials, |J3], one then observes that the orbits do indeed disperse in such 
a fashion as to exhibit a coarse-grained evolution towards a quasi-equilibrium, which 
is at least approximately time-independent. This approach is, moreover, exponential 
in time and characterised by a rate A(E) which is independent of the specific choice 
of initial data. The characteristic timescale t* = A -1 is typically <C 100t cr , so that, 
in "physical units" for a galaxy, <C 

One also observes that the rate A(E) is comparable in magnitude to the value 
of the Liapounov exponent x{E), which|35| probes the average rate of instability 
exhibited by orbits of energy E. There is, moreover, a direct correlation between A 
and x m the sense that, e.g., both have similar curvatures when viewed as functions of 
E. This is particularly tantalising in view of the fact that, for the A^-body problem, 
the timescale associated with the approach towards a statistical quasi-equilibrium 
is comparable in magnitude to the timescale associated with the instability towards 
small changes in initial conditions. 

Despite these regularities, there is no guarantee that this apparent equilibrium 
coincides with the true invariant distribution, and, in general, it will not\ Astronomers 
are well acquainted with the fact that collisionless equilibria do not constitute true N- 
body equilibria, which must incorporate discreteness effects that become important 
on sufficiently long timescales. However, the key point here is different, and more 
fundamental: Even motion in a smooth two-dimensional potential may not yield a 
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uniform approach towards a true statistical equilibrium, ,[34. 36 

The explanation of the discrepancy between the true and approximate equilibria is 
simple. Viewed over sufficiently long timescales, there are only two different classes of 
orbits, namely regular orbits, with vanishing Liapounov exponent x> an d stochastic 
orbits, for which x > 0. The distinction between these two classes is, moreover, 
absolute, since members of the two different classes are separated by invariant KAM 
tori. If, e.g., one were to compute a surface of section, plotting coordinates x and p x 
for successive intersections of the y = hyperplane, he or she would generically find 
islands of regularity embedded in a surrounding stochastic sea. 

However, this is not the whole story. Lurking in the shallows of the stochastic 
sea, slightly away from the shore, are cantori, |37]] these corresponding to fractured 
KAM tori, associated with the breakdown of integrability, which contain a cantor set 
of holes. The point then is that these cantori serve as partial barriers that divide 
the stochastic orbits into two subclasses, namely confined, or sticky, stochastic orbits 
which are trapped near the regular islands, and unconfined, or filling, stochastic orbits 
which travel unimpeded throughout the rest of the stochastic sea. 

Because of the holes in the cantori, these barriers are not absolute, so that orbits 
can in fact change from one class to another via so-called intrinsic diffusion. However, 
this process is a slow one, requiring orbits to wend their way through a maze (cf. 



the "turnstile model" of MacKay, Meiss, and Percival,[^] so that the characteristic 
timescale is typically ^> 100t cr , i.e., much longer than the age of the Universe. 

What this implies is that, on short times, ensembles of orbits initially outside the 
cantori will evolve towards a near-invariant distribution which uniformly populates 
the filling regions, but avoids the confined regions. The situation is analogous to 
the classical effusion problem. Consider two evacuated cavities connected one with 
another by an extremely narrow conduit, and suppose that gas is inserted into one of 
the cavities. If the conduit be sufficiently narrow, the timescale on which gas effuses 
from one cavity to the other will be much longer than the timescale on which the gas 
spreads to fill the original cavity. This implies, however, that, even though the true 
equilibrium corresponds to a uniform density concentration throughout both cavities, 
one can speak meaningful of a shorter time quasi-equilibrium, in which the original 
cavity is populated uniformly and the other is essentially empty. 

Significantly, these two different populations of stochastic orbits are fundamentally 
dissimilar in terms of their stability properties as well as where they are located in 
phase space. Although both sticky and filling stochastic orbits are exponentially 
unstable, there is a precise sense in which the sticky orbits are less unstable overall 



than are the filling orbits. Specifically, if one computes local Liapounov exponents, \ 39 1 
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x(At), for different ensembles of stochastic orbits, integrating for some relatively short 
interval At, he or she will findp4|, |36| that the typical x(At) for a sticky orbit is 
substantially smaller than the typical xi^t) for a filling orbit. Indeed, the composite 
distribution of local Liapounov exponents (i.e., distribution of instability timescales) 
generated from a sampling of the true invariant measure appears to be given, at least 
approximately, as a sum of two different near-Gaussian distributions with unequal 
means. 

It should perhaps be noted explicitly that the general conclusions recounted in 
this Section have been observed for several different potentials, with rather different 
symmetries, including (1) the dihedral D4 potential of Armbruster, Guckenheimer, 
and Kim,^TJ (2) the sixth order truncation of the three-particle Toda|4l| lattice 
potential, and (3) a generalised anisotropic Kepler potential of the form 

1 TYl 

V ( x ' V) = ~7 ZTJ2 - 7 7T72> ( 7 ) 

(l + x 2 + y 2 ) (l + x 2 + ay 2 ) 

with constant m and a, for E < 0. The fact that these diverse potentials, which are 
fundamentally different in appearance, yield similar conclusions, both qualitatively 
and semi-quantitatively, would suggest strongly that these conclusions are robust, 
depending only on such topological features as the existence of KAM tori and cantori. 

The existence of confined stochastic orbits is of potential importance astronomi- 
cally because such orbits can help (the theorist) support various sorts of structures, 
e.g., bars in a spiral galaxy. It is natural to assume that, in systems like galaxies, 
regular orbits serve to provide the skeleton to support various structures. However, 
because of resonance overlap one may find that, near corotation and other resonances, 
the desired regular orbits do not exist, even though sticky stochastic orbits are present. 

Finally, it should be stressed that one can observe similar short time "zones of 
avoidance" in higher dimensional systems as well. The key point physically is that 
just because a region of phase space is connected, so that orbits can pass throughout 
the entire region, does not mean that all of the region will be accessed on comparable 
timescales. 



5 Structural stability of the smooth potential ap- 
proximation 

The collisionless Boltzmann equation is a Hamiltonian system which neglects various 
realistic non-Hamiltonian irregularities that must be present in any self-gravitating 
system. One obvious point is that such a Vlasov description neglects entirely all 
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discreteness effects, i.e., "collisions," by idealising the system as a continuum, rather 
than a collection of nearly point mass objects. Viewed in the iV-particle phase space, 
the statistical description of an isolated iV-body evolution is of course Hamiltonian. 
However, when projected into the reduced one-particle phase space, any allowance 
for particle-particle correlations that transcend a mean field description necessar- 
ily breaks the Hamiltonian constraints.^] Another point, perhaps less obvious but 
equally important, is that a Vlasov description also neglects any couplings to an ex- 
ternal environment. In the past, astronomers have been wont oftentimes to pretend 
that galaxies exist in splendid isolation but, over the past several decades, it has 
become increasing evident that such an approximation may not be justified. [fl2f 

Detailed modeling of these sorts of perturbing influences may prove extremely 
complex. In particular, an external environment can give rise to a variety of different 
effects characterised by a broad range of timescales. Those influences proceeding on 
timescales ~ t cr will be particularly complicated, in that the details of their effects 
may depend very sensitively on the details of the environment. However, there is 
a well-established paradigm in statistical physics, f43|, |44}| dating back to the begin- 



ning of the century, ||45|| which would suggest that irregularities proceeding on shorter 



timescales, <C t cr , can oftentimes be modelled as friction and noise, related via a 
fluctuation-dissipation theorem. This idea underlies, for example, Chandrasekhar's 
original formulation [f5| of so-called "collisional stellar dynamics." 



It is therefore natural to investigate the structural stability of Hamiltonian tra- 



jectories towards the effects of friction and noise. This was done [f47| , [48] by effecting 
large numbers of Langevin simulations, in which the deterministic equations of motion 
were perturbed by allowing for (1) a dynamical friction —f]p, which serves system- 
atically to remove energy from the orbits and (2) random kicks, modeled as white 
noise with temperature, or mean squared velocity, 0, which serve systematically to 
pump energy back into the orbits. As a first simple test, r] was assumed to be con- 
stant, in which case the fluctuation-dissipation theorem implies that the noise must 
be additive, rather than multiplicative.^, ^TJ 

Thus, in units with particle mass m = 1, one is led explicitly to equations of 
motion of the form 

Tt =P * = -Vy(r)- W + F, (8) 

where F is characterised completely by its statistical properties. Here, e.g., component 
by component, 



(Fi(t)) = and (F^Fjih)) = 2Qr ] S ij S D (t 1 - t 2 ), 
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(9) 



where the angular brackets denote a statistical average. The idea is to effect large 
numbers of different realisations of the same initial conditions, and to analyse these 
realisations to extract statistical properties. 

It is well known that even very weak friction and noise will eventually become 
important on sufficiently long timescales. In particular, one knows that, on the natural 
timescale ~ r]~ l , these effects will try to force the system to evolve towards a 
thermal state. The question of relevance here is quite different: Can the friction 
and noise have substantial effects already on much shorter timescales -C £r? The 
conventional wisdom in astronomy is that the answer to this is: no\ For example, the 
standard assumption that "collisionality" is irrelevant in a galaxy relies completely 
on the observation that the natural timescale associated with binary encounters is 



much longer than the age of the Universe. |4"6[| 



The Langevin simulations were effected for total times t < 200t cr , and involved 
friction and noise corresponding to a broad range of characteristic timescales, 10 3 < 
tR/t cr < 10 12 . The most significant conclusions derived from these simulations are 
the following: 

When viewed in terms of the collisionless invariants, i.e., the quantities that are 
conserved in the absence of the friction and noise, these perturbing influences only 
serve to induce a classical diffusion process, with the unperturbed and perturbed 
orbits diverging significantly only on a timescale t# ~ r/ _1 . Thus, in particular, 

SE 2 rms = (\E unp - E per \ 2 ) = A(E)EQrjt, (10) 

where A(E) is a slowly varying function of E with magnitude of order unity. In this 
sense, the conventional wisdom is confirmed. 

However, when viewed in configuration space or momentum space, the effects are 
more complicated, and actually depend on orbit class. Unperturbed and perturbed 
regular orbits only diverge as a power law in time, i.e., 5r rms , 6p rms ~ t p , so that, 
once again, one only gets macroscopic deviations after a time tu ~ r/^ 1 . However, 
unperturbed and perturbed stochastic orbits diverge exponentially at a rate set by 
the Liapounov exponent Xi so that, even for very weak friction and noise, one gets 
macroscopic deviations within a few crossing times. In particular, when considering 
ensembles of stochastic initial conditions, one observes a simple scaling 

5r rms , 5 Prms ~ (e??) 1/2 exp[+X(£)t], (11) 

where X is comparable to, but slightly larger than, the Liapounov exponent x- 

This exponential divergence is easy to understand p2[ and the functional depen- 



dence on 9, rj, and x can actually be derived theoretically. j53| The obvious point is 
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that the unperturbed deterministic trajectory is an unstable stochastic orbit, so that 
even the tiniest perturbing influences will tend to grow exponentially. The average 
rate of instability is given by x an d, as such, moments like Sr rms should grow at a rate 
X ~ x- That X is slightly larger than x is a reflection of the fact that, for different 
noisy realisations, one sees somewhat different local Liapounov exponents, and that 
the total Sr rms will be dominated by those noisy realisations for which the rate of 
instability is above average. 

This argument might suggest that, although the unperturbed and perturbed tra- 
jectories exhibit a rapid pointwise divergence, their statistical properties should be 
virtually identical. Specifically, one might anticipate that, on short times, the only ef- 
fect of the friction and noise is to continually displace the trajectory from one stochas- 
tic orbit to another with essentially the same statistical properties. This, however, is 
false. Under certain circumstances, even very weak friction and noise can also alter 
the statistical properties of ensembles of stochastic orbits on relatively short times 
<C 100t cr . Specifically, one observes that such perturbing influences can dramatically 
accelerate the rate of penetration through cantori by providing an additional source 
of extrinsic diffusion. 

Provided that the friction and noise are sufficiently weak, on short timescales the 
energy E is almost conserved, so that one can speak meaningfully of an evolution 
restricted to an "almost constant energy hypersurface." Suppose now that, for some 
energy E, this hypersurface contains large measures of both sticky and filling stochas- 
tic orbits. If, for this energy, the near-invariant distribution described in Section 4 
is evolved into the future, allowing for even very weak friction and noise, one then 
observes a rapid (t <C 100t cr ) systematic evolution towards a new noisy near- invariant 
distribution which is (1) quite different from the deterministic near-invariant distribu- 
tion and (2) much closer to the true deterministic invariant distribution. In this sense, 
it appears that, on timescales short compared the timescale t^ on which the system 
would evolve towards a thermal state, the principal effect of the friction and noise is 
to accelerate the approach towards a deterministic invariant distribution which, in the 
absence of these perturbing influences, would only have been realised on much longer 
timescales. 

The key point in all of this is that friction and noise can induce changes in orbit 
class, from filling to confined stochastic, and vice versa. Moreover, when the deter- 
ministic invariant distribution contains large measures of both sticky and filling orbits, 
such transitions can happen within a time t < 100t cr , even for very weak friction and 
noise. Visual inspection of ~ 2.5 x 10 4 orbits in several different potentials leads to 
the following conclusions. 
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Typically, for a relaxation time tn as long as 10 12 t cr , not many such changes are 
observed within a time t ~ 100t cr . However, if tn be reduced to a value ~ 10 9 t cr , 
transitions begin to become more frequent, and, even for t^ as large as ~ 10 6 t 
transitions are quite common, occuring for > 50% of orbits within a time t ~ lOOt 
If the amplitude of the friction and noise are further increased, one finds that, for 

~ 10 3 t cr , transitions are so common that the distinction between filling and con- 
fined becomes essentially meaningless. The distinction between regular and stochastic 
is more robust. Only for Ir as small as ~ 10 3 t cr are significant numbers of transitions 
between regular and stochastic orbits observed within a time as short as t ~ 100t cr . 

The fact that even very weak friction and noise, with tjt ~ 10 6 t cr — 10 9 t cr , can 
significantly alter the statistical properties of ensembles of orbits on timescales t < 
tu ~ 100t cr has direct astronomical implications since, e.g., for galaxies, the timescale 
in for discreteness effects, i.e., "collisionality" is ~ 10 6 t cr ! The natural timescale 
associated with external perturbations is less easily estimated, but may well be even 
shorter. 

To summarise, it is evident that even very weak friction and noise can alter both 
the pointwise and the statistical properties of stochastic orbits in a nonintegrable 
potential on relative short timescales C^. In particular, such effects may be man- 
ifest already on time scales much shorter than the time on which numerical errors 
in a simulation can accumulate. This fact has direct and immediate implications 



for the problem of "shadowing" for numerical orbits. |54|, |55| Physicists, mathemati- 



cians, and astronomers are often worried [56], 57 1 about whether numerical simulations 
performed on a computer, which incorporate roundoff and/or truncation error, can 
correctly shadow the evolution of some model system described by a simple set of de- 
terministic differential equations. However, it would also seem relevant |58| to worry 
about whether the "real world," replete with other sorts of irregularities, can shadow 
either the model system or its numerical realisations. In this regard, one final remark 
is in order: Rather than being an impediment to realistic modeling, in certain cases 
numerical noise may actually be a good thing, in that it may capture, at least quali- 
tatively, some of the effects of small perturbing influences to which real systems are 
always subjected. 
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